COVID19 UK

Prepare data

Read and format data

Prevalence

df_uk_covid$time %>% summary()

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.00    8.75   16.50   16.50   24.25   32.00 

Social distancing

df_uk_socdist$date %>% summary()

        Min.      1st Qu.       Median         Mean      3rd Qu.         Max. 
"2020-03-01" "2020-03-08" "2020-03-16" "2020-03-16" "2020-03-24" "2020-03-31" 

Controls

df_uk_ctrl_nuts <- read_csv("controls_UK_nuts3.csv")

Parsed with column specification:
cols(
  nuts3 = col_character(),
  nuts3_name = col_character(),
  airport_dist = col_double(),
  males = col_double(),
  popdens = col_double(),
  manufacturing = col_double(),
  tourism = col_double(),
  health = col_double(),
  academic = col_double(),
  medinc = col_double(),
  medage = col_double(),
  conservative = col_double()
)

df_uk_ctrl_nuts <- df_uk_ctrl_nuts %>% select(-nuts3_name)
df_uk_ctrl_nuts

df_uk_ctrl_ut <- read_csv("controls_UK_ut.csv")

Parsed with column specification:
cols(
  ut_area = col_character(),
  ut_name = col_character(),
  airport_dist = col_double(),
  males = col_double(),
  popdens = col_double(),
  manufacturing = col_double(),
  tourism = col_double(),
  health = col_double(),
  academic = col_double(),
  medinc = col_double(),
  medage = col_double(),
  conservative = col_double()
)

df_uk_ctrl_ut <- df_uk_ctrl_ut %>% select(-ut_name)
df_uk_ctrl_ut

NA
NA

Merge data

df_uk <- df_uk %>% plyr::join(df_uk_ctrl_ut, by='ut_area')

df_uk_socdist <- df_uk_socdist %>% plyr::join(df_uk_ctrl_nuts, by='nuts3')

Identify London areas

nuts_london_inner <- c('UKI31','UKI32','UKI33','UKI34','UKI41',
                      'UKI42','UKI43','UKI44','UKI45')

nuts_london_outer <- c('UKI51','UKI52','UKI53','UKI54','UKI61',
                      'UKI62','UKI63','UKI71','UKI72','UKI73',
                      'UKI74','UKI75')

ut_london_inner <- c('E09000007','E09000001','E09000033','E09000013',
                    'E09000020','E09000032','E09000025','E09000012',
                    'E09000030','E09000014','E09000019','E09000023',
                    'E09000028','E09000022')

ut_london_outer <- c('E09000011','E09000004','E09000016','E09000002',
                    'E09000031','E09000026','E09000010','E09000006',
                    'E09000008','E09000029','E09000021','E09000024',
                    'E09000003','E09000005','E09000009','E09000017',
                    'E09000015','E09000018','E09000027')

df_uk = df_uk %>% 
  mutate(london = ifelse(ut_area %in% ut_london_inner, 'london_inner', 
                       ifelse(ut_area %in% ut_london_outer, 'london_outer',
                              'country'))) %>%
  mutate(london = as.factor(london))

df_uk_socdist = df_uk_socdist %>% 
  mutate(london = ifelse(nuts3 %in% nuts_london_inner, 'london_inner', 
                       ifelse(nuts3 %in% nuts_london_outer, 'london_outer',
                              'country'))) %>%
  mutate(london = as.factor(london))

Explore data

Plot prevalence over time

df_uk %>% ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=ut_area, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall prevalence over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_uk %>% mutate(prev_tail = cut(.[[i]], 
                                       breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
                                       labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(prev_tail != 'center') %>%
  ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=ut_area, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~prev_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

Plot social distancing over time

df_uk_socdist %>% ggplot(aes(x=time, y=socdist_single_tile)) + 
  geom_point(aes(col=nuts3, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall social distancing over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_uk_socdist %>% mutate(socdist_tail = cut(.[[i]], 
                                       breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
                                       labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(socdist_tail != 'center') %>%
  ggplot(aes(x=time, y=socdist_single_tile)) + 
  geom_point(aes(col=nuts3, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~socdist_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

Explore differences between london and the rest

df_uk %>% ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=ut_area, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  facet_wrap(~london) +
  ggtitle("Overall prevalence over time")

NA
NA
NA

df_uk_socdist %>% ggplot(aes(x=time, y=socdist_single_tile)) + 
  geom_point(aes(col=nuts3, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  facet_wrap(~london) +
  ggtitle("Overall social distancing over time")

Control for weekend effect

weekend <- c(6, 7, 13, 14, 20, 21)

df_uk_loess <- df_uk_socdist %>% filter(!time %in% weekend) %>% 
  split(.$nuts3) %>%
  map(~ loess(socdist_single_tile ~ time, data = .)) %>%
  map(predict, 1:23) %>% 
  bind_rows() %>% 
  gather(key = 'nuts3', value = 'loess') %>% 
  group_by(nuts3) %>% 
  mutate(time = row_number())

df_uk_socdist <- df_uk_socdist %>% merge(df_uk_loess, by=c('nuts3', 'time')) %>% 
  mutate(socdist_single_tile_clean = ifelse(time %in% weekend, loess,
                                            socdist_single_tile)) %>%
  arrange(nuts3, time)


df_uk_socdist %>% ggplot(aes(x=time, y=loess, group=nuts3)) +
  geom_line()

df_uk_socdist

NA
NA

df_uk_socdist %>% ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=nuts3, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall social distancing over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_uk_socdist %>% mutate(socdist_tail = cut(.[[i]], 
                                       breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
                                       labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(socdist_tail != 'center') %>%
  ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=nuts3, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~socdist_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

df_uk_socdist <- df_uk_socdist %>% mutate(socdist_single_tile = socdist_single_tile_clean) %>% 
  select(-loess, -socdist_single_tile_clean)

Correlations

df_uk %>% group_by(ut_area) %>% 
  summarize_if(is.numeric, mean, na.rm=T) %>% 
  select(-ut_area, -time) %>% 
  cor(use = 'pairwise.complete') %>% round(3)

              pers_o pers_e pers_a pers_n pers_c  frequ rate_day airport_dist  males popdens
pers_o         1.000  0.734 -0.604 -0.143 -0.636  0.011    0.647       -0.247  0.437   0.806
pers_e         0.734  1.000 -0.412 -0.497 -0.293  0.071    0.559       -0.262  0.197   0.617
pers_a        -0.604 -0.412  1.000 -0.197  0.643  0.130   -0.584        0.285 -0.460  -0.745
pers_n        -0.143 -0.497 -0.197  1.000 -0.353 -0.243   -0.108        0.059  0.096   0.043
pers_c        -0.636 -0.293  0.643 -0.353  1.000  0.220   -0.487        0.327 -0.524  -0.727
frequ          0.011  0.071  0.130 -0.243  0.220  1.000   -0.077        0.108 -0.124  -0.180
rate_day       0.647  0.559 -0.584 -0.108 -0.487 -0.077    1.000       -0.388  0.302   0.758
airport_dist  -0.247 -0.262  0.285  0.059  0.327  0.108   -0.388        1.000 -0.124  -0.381
males          0.437  0.197 -0.460  0.096 -0.524 -0.124    0.302       -0.124  1.000   0.510
popdens        0.806  0.617 -0.745  0.043 -0.727 -0.180    0.758       -0.381  0.510   1.000
manufacturing -0.532 -0.595  0.496  0.359  0.363 -0.111   -0.511        0.300 -0.259  -0.545
tourism        0.116  0.098  0.069 -0.244  0.227  0.115   -0.023        0.448 -0.094  -0.099
health         0.089  0.026 -0.001  0.138 -0.286 -0.132   -0.010        0.150 -0.012   0.112
academic       0.707  0.730 -0.502 -0.447 -0.298  0.152    0.647       -0.369  0.267   0.602
medinc         0.533  0.569 -0.533 -0.335 -0.300  0.070    0.573       -0.379  0.367   0.631
medage        -0.508 -0.337  0.582 -0.167  0.780  0.144   -0.513        0.447 -0.620  -0.705
conservative  -0.843 -0.767  0.511  0.366  0.485 -0.125   -0.626        0.396 -0.293  -0.686
              manufacturing tourism health academic medinc medage conservative
pers_o               -0.532   0.116  0.089    0.707  0.533 -0.508       -0.843
pers_e               -0.595   0.098  0.026    0.730  0.569 -0.337       -0.767
pers_a                0.496   0.069 -0.001   -0.502 -0.533  0.582        0.511
pers_n                0.359  -0.244  0.138   -0.447 -0.335 -0.167        0.366
pers_c                0.363   0.227 -0.286   -0.298 -0.300  0.780        0.485
frequ                -0.111   0.115 -0.132    0.152  0.070  0.144       -0.125
rate_day             -0.511  -0.023 -0.010    0.647  0.573 -0.513       -0.626
airport_dist          0.300   0.448  0.150   -0.369 -0.379  0.447        0.396
males                -0.259  -0.094 -0.012    0.267  0.367 -0.620       -0.293
popdens              -0.545  -0.099  0.112    0.602  0.631 -0.705       -0.686
manufacturing         1.000  -0.088 -0.023   -0.714 -0.532  0.472        0.657
tourism              -0.088   1.000  0.049    0.056 -0.062  0.435       -0.023
health               -0.023   0.049  1.000   -0.228 -0.243 -0.153        0.043
academic             -0.714   0.056 -0.228    1.000  0.733 -0.382       -0.888
medinc               -0.532  -0.062 -0.243    0.733  1.000 -0.481       -0.661
medage                0.472   0.435 -0.153   -0.382 -0.481  1.000        0.486
conservative          0.657  -0.023  0.043   -0.888 -0.661  0.486        1.000

df_uk_socdist %>% group_by(nuts3) %>% 
  summarize_if(is.numeric, mean, na.rm=T) %>% 
  select(-nuts3, -time) %>% 
  cor(use = 'pairwise.complete') %>% round(3)

                    rate_day socdist_tiles socdist_single_tile pers_o pers_e pers_a pers_n pers_c  frequ
rate_day               1.000        -0.695               0.487  0.635  0.562 -0.567 -0.113 -0.469  0.078
socdist_tiles         -0.695         1.000              -0.619 -0.626 -0.733  0.545  0.325  0.395 -0.311
socdist_single_tile    0.487        -0.619               1.000  0.317  0.356 -0.382  0.040 -0.298  0.032
pers_o                 0.635        -0.626               0.317  1.000  0.713 -0.603 -0.122 -0.651  0.183
pers_e                 0.562        -0.733               0.356  0.713  1.000 -0.450 -0.456 -0.363  0.249
pers_a                -0.567         0.545              -0.382 -0.603 -0.450  1.000 -0.215  0.661 -0.030
pers_n                -0.113         0.325               0.040 -0.122 -0.456 -0.215  1.000 -0.338 -0.332
pers_c                -0.469         0.395              -0.298 -0.651 -0.363  0.661 -0.338  1.000  0.096
frequ                  0.078        -0.311               0.032  0.183  0.249 -0.030 -0.332  0.096  1.000
airport_dist          -0.353         0.550              -0.357 -0.231 -0.319  0.328  0.072  0.326 -0.147
males                  0.303        -0.287               0.077  0.519  0.243 -0.567  0.195 -0.632 -0.032
popdens                0.730        -0.705               0.500  0.803  0.638 -0.748  0.030 -0.737 -0.040
manufacturing         -0.412         0.668              -0.438 -0.496 -0.605  0.427  0.343  0.309 -0.278
tourism               -0.048         0.239              -0.110  0.079  0.005  0.142 -0.168  0.207 -0.007
health                 0.002         0.140              -0.060  0.086  0.047  0.005  0.151 -0.274 -0.258
academic               0.640        -0.808               0.348  0.713  0.742 -0.491 -0.432 -0.324  0.427
medinc                 0.600        -0.772               0.457  0.552  0.582 -0.620 -0.247 -0.367  0.248
medage                -0.472         0.563              -0.328 -0.512 -0.402  0.631 -0.181  0.790 -0.047
conservative          -0.629         0.757              -0.327 -0.837 -0.761  0.514  0.338  0.522 -0.374
                    airport_dist  males popdens manufacturing tourism health academic medinc medage
rate_day                  -0.353  0.303   0.730        -0.412  -0.048  0.002    0.640  0.600 -0.472
socdist_tiles              0.550 -0.287  -0.705         0.668   0.239  0.140   -0.808 -0.772  0.563
socdist_single_tile       -0.357  0.077   0.500        -0.438  -0.110 -0.060    0.348  0.457 -0.328
pers_o                    -0.231  0.519   0.803        -0.496   0.079  0.086    0.713  0.552 -0.512
pers_e                    -0.319  0.243   0.638        -0.605   0.005  0.047    0.742  0.582 -0.402
pers_a                     0.328 -0.567  -0.748         0.427   0.142  0.005   -0.491 -0.620  0.631
pers_n                     0.072  0.195   0.030         0.343  -0.168  0.151   -0.432 -0.247 -0.181
pers_c                     0.326 -0.632  -0.737         0.309   0.207 -0.274   -0.324 -0.367  0.790
frequ                     -0.147 -0.032  -0.040        -0.278  -0.007 -0.258    0.427  0.248 -0.047
airport_dist               1.000 -0.155  -0.377         0.272   0.492  0.180   -0.379 -0.389  0.469
males                     -0.155  1.000   0.615        -0.180  -0.171  0.028    0.263  0.440 -0.672
popdens                   -0.377  0.615   1.000        -0.492  -0.135  0.168    0.574  0.652 -0.714
manufacturing              0.272 -0.180  -0.492         1.000  -0.051 -0.069   -0.640 -0.462  0.413
tourism                    0.492 -0.171  -0.135        -0.051   1.000  0.027    0.005 -0.138  0.485
health                     0.180  0.028   0.168        -0.069   0.027  1.000   -0.202 -0.261 -0.154
academic                  -0.379  0.263   0.574        -0.640   0.005 -0.202    1.000  0.724 -0.368
medinc                    -0.389  0.440   0.652        -0.462  -0.138 -0.261    0.724  1.000 -0.479
medage                     0.469 -0.672  -0.714         0.413   0.485 -0.154   -0.368 -0.479  1.000
conservative               0.389 -0.350  -0.678         0.627   0.012  0.010   -0.891 -0.653  0.499
                    conservative
rate_day                  -0.629
socdist_tiles              0.757
socdist_single_tile       -0.327
pers_o                    -0.837
pers_e                    -0.761
pers_a                     0.514
pers_n                     0.338
pers_c                     0.522
frequ                     -0.374
airport_dist               0.389
males                     -0.350
popdens                   -0.678
manufacturing              0.627
tourism                    0.012
health                     0.010
academic                  -0.891
medinc                    -0.653
medage                     0.499
conservative               1.000

Modelling

Prepare functions

# function calculates all relevant models
run_models <- function(y, lvl1_x, lvl2_x, lvl2_id, data, ctrls=F){

  # subset data
  data = data %>% 
    dplyr::select(all_of(y), all_of(lvl1_x), all_of(lvl2_x), all_of(lvl2_id), 
                  popdens, rate_day, all_of(y))
  data = data %>% 
    dplyr::rename(y = all_of(y),
           lvl1_x = all_of(lvl1_x),
           lvl2_x = all_of(lvl2_x),
           lvl2_id = all_of(lvl2_id)
           )
  
  # configure optimization procedure
  ctrl_config <- lmeControl(opt = 'optim', maxIter = 100, msMaxIter = 100)

  # baseline
  baseline <- lme(fixed = y ~ 1, random = ~ 1 | lvl2_id, 
                    data = data,
                    correlation = corAR1(),
                  control = ctrl_config,
                  method = 'ML')

  # random intercept fixed slope
  random_intercept <- lme(fixed = y ~ lvl1_x + lvl2_x, 
                          random = ~ 1 | lvl2_id,
                            data = data,
                            correlation = corAR1(),
                  control = ctrl_config,
                  method = 'ML')

  # random intercept random slope
  random_slope <- lme(fixed = y ~ lvl1_x + lvl2_x, 
                      random = ~ lvl1_x | lvl2_id, 
                        data = data,
                        correlation = corAR1(),
                  control = ctrl_config,
                  method = 'ML')

  # cross level interaction
  interaction <- lme(fixed = y ~ lvl1_x * lvl2_x, 
                     random = ~ lvl1_x | lvl2_id, 
                       data = data,
                       correlation = corAR1(),
                  control = ctrl_config,
                  method = 'ML')
  
  # create list with results
  results <- list('baseline' = baseline, 
                  "random_intercept" = random_intercept, 
                  "random_slope" = random_slope,
                  "interaction" = interaction)
  
  
  if (ctrls == 'dem' | ctrls == 'prev'){
    
    # random intercept random slope
    random_slope_ctrl_dem <- lme(fixed = y ~ lvl1_x + lvl2_x + popdens,
                              random = ~ lvl1_x | lvl2_id, 
                          data = data,
                          correlation = corAR1(),
                    control = ctrl_config,
                  method = 'ML')
  
    # cross level interaction
    interaction_ctrl_main_dem <- lme(fixed = y ~ lvl1_x * lvl2_x + popdens,
                             random = ~ lvl1_x | lvl2_id, 
                         data = data,
                         correlation = corAR1(),
                    control = ctrl_config,
                  method = 'ML')
  
    # cross level interaction
    interaction_ctrl_int_dem <- lme(fixed = y ~ lvl1_x * lvl2_x + lvl1_x * popdens,
                             random = ~ lvl1_x | lvl2_id, 
                         data = data,
                         correlation = corAR1(),
                    control = ctrl_config,
                  method = 'ML')        
    
    # create list with results
    results <- list('baseline' = baseline, 
                    "random_intercept" = random_intercept, 
                    "random_slope" = random_slope,
                    "interaction" = interaction,
                    "random_slope_ctrl_dem" = random_slope_ctrl_dem,
                    "interaction_ctrl_main_dem" = interaction_ctrl_main_dem,
                    "interaction_ctrl_int_dem" = interaction_ctrl_int_dem)
  }
  
  if (ctrls == 'prev'){
  
    # random intercept random slope
    random_slope_ctrl_prev <- lme(fixed = y ~ lvl1_x + lvl2_x + popdens + rate_day,
                              random = ~ lvl1_x + rate_day | lvl2_id, 
                          data = data,
                          correlation = corAR1(),
                          control = ctrl_config,
                  method = 'ML')  
    
        # cross level interaction
    interaction_ctrl_main_prev <- lme(fixed = y ~ lvl1_x * lvl2_x + popdens + rate_day,
                             random = ~ lvl1_x | lvl2_id, 
                         data = data,
                         correlation = corAR1(),
                    control = ctrl_config,
                  method = 'ML')
  
  
    # cross level interaction
    interaction_ctrl_int_prev<- lme(fixed = y ~ lvl1_x * lvl2_x + lvl1_x * popdens + rate_day,
                             random = ~ lvl1_x + rate_day | lvl2_id, 
                         data = data,
                         correlation = corAR1(),
                          control = ctrl_config,
                  method = 'ML')
  
    # create list with results
    results <- list('baseline' = baseline, 
                    "random_intercept" = random_intercept, 
                    "random_slope" = random_slope,
                    "interaction" = interaction,
                    "random_slope_ctrl_dem" = random_slope_ctrl_dem,
                    "interaction_ctrl_main_dem" = interaction_ctrl_main_dem,
                    "interaction_ctrl_int_dem" = interaction_ctrl_int_dem,                    
                    "random_slope_ctrl_prev" = random_slope_ctrl_prev,
                    "interaction_ctrl_main_prev" = interaction_ctrl_main_prev,
                    "interaction_ctrl_int_prev" = interaction_ctrl_int_prev)
  }
  
  if(ctrls == 'exp'){
    # random intercept random slope
  random_slope_exp <- lme(fixed = y ~ (lvl1_x + I(lvl1_x^2)) + lvl2_x, 
                      random = ~ (lvl1_x + I(lvl1_x^2)) | lvl2_id, 
                        data = data,
                        correlation = corAR1(),
                  control = ctrl_config,
                  method = 'ML')

  # cross level interaction
  interaction_exp <- lme(fixed = y ~ (lvl1_x + I(lvl1_x^2)) * lvl2_x, 
                     random = ~ (lvl1_x + I(lvl1_x^2)) | lvl2_id, 
                       data = data,
                       correlation = corAR1(),
                  control = ctrl_config,
                  method = 'ML')  
  
  
  # create list with results
  results <- list('baseline' = baseline, 
                  "random_intercept" = random_intercept, 
                  "random_slope" = random_slope,
                  "interaction" = interaction,                  
                  "random_slope_exp" = random_slope_exp,
                  "interaction_exp" = interaction_exp)
  }
  
  return(results)
        
}

# extracts table with coefficients and tests statistics
extract_results <- function(models) {
  
  models_summary <- models %>% 
  map(summary) %>% 
  map("tTable") %>% 
  map(as.data.frame) %>% 
  map(round, 10) 
  # %>% map(~ .[str_detect(rownames(.), 'Inter|lvl'),])
  
  return(models_summary)
  
}


# calculates comparison of all models in model list
compare_models <- function(models) {

  mdl_names <- models %>% names()
  
  str = ''
  for (i in mdl_names){
    
    mdl_str <- paste('models$', i, sep = '')
    
    if(i == 'baseline'){
      str <- mdl_str
    }else{
    str <- paste(str, mdl_str, sep=', ')
    }
  }
  
  anova_str <- paste0('anova(', str, ')')
  mdl_comp <- eval(parse(text=anova_str))
  rownames(mdl_comp) = mdl_names
  return(mdl_comp)
}

Remove London Data

# df_uk <- df_uk %>% filter(london == 'country')
# df_uk_socdist <- df_uk_socdist %>% filter(london == 'country')

Rescale Data

lvl2_scaled_ut <- df_uk %>% 
  dplyr::select(-time, -frequ, -rate_day, -london) %>% 
  distinct() %>% 
  mutate_at(vars(-ut_area), scale)

lvl1_scaled_ut <- df_uk %>% select(ut_area, time, rate_day) %>% 
  mutate_at(vars(-ut_area, -time), scale)

df_uk_scaled <- plyr::join(lvl1_scaled_ut, lvl2_scaled_ut, by = 'ut_area')

head(df_uk_scaled)

lvl2_scaled_nuts <- df_uk_socdist %>% 
  dplyr::select(-time, -date, -frequ, -london, 
                -socdist_tiles, -socdist_single_tile, -rate_day) %>% 
  distinct() %>% 
  mutate_at(vars(-nuts3), scale)

lvl1_scaled_nuts <- df_uk_socdist %>% select(nuts3, time, socdist_single_tile, rate_day) %>% 
  mutate_at(vars(-nuts3, -time), scale)

df_uk_socdist_scaled <- plyr::join(lvl1_scaled_nuts, lvl2_scaled_nuts, by = 'nuts3')

head(df_uk_socdist_scaled)

NA

Predict prevalence

prevalence ~ openness

models_o_covid <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_o', 
                         lvl2_id = 'ut_area', 
                         data = df_uk_scaled,
                         ctrls = 'dem')

extract_results(models_o_covid)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

compare_models(models_o_covid)

NA

prevalence ~ conscientiousness

models_c_covid <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_c', 
                         lvl2_id = 'ut_area', 
                         data = df_uk_scaled,
                         ctrls = 'dem')

extract_results(models_c_covid)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

compare_models(models_c_covid)

NA
NA

prevalence ~ extraversion

models_e_covid <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_e', 
                         lvl2_id = 'ut_area', 
                         data = df_uk_scaled,
                         ctrls = 'dem')

extract_results(models_e_covid)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

compare_models(models_e_covid)

NA
NA

prevalence ~ agreeableness

models_a_covid <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_a', 
                         lvl2_id = 'ut_area', 
                         data = df_uk_scaled,
                         ctrls = 'dem')

extract_results(models_a_covid)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

compare_models(models_a_covid)

NA
NA

prevalence ~ neuroticism

models_n_covid <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_n', 
                         lvl2_id = 'ut_area', 
                         data = df_uk_scaled,
                         ctrls = 'dem')

extract_results(models_n_covid)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

compare_models(models_n_covid)

NA
NA

Predict social distancing

social distancing ~ openness

models_o_socdist <-run_models(y = 'socdist_single_tile', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_o', 
                         lvl2_id = 'nuts3', 
                         data = df_uk_socdist_scaled,
                         ctrls = 'prev')

extract_results(models_o_socdist)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

$random_slope_ctrl_prev

$interaction_ctrl_main_prev

$interaction_ctrl_int_prev

compare_models(models_o_socdist)

NA

social distancing ~ conscientiousness

models_c_socdist <-run_models(y = 'socdist_single_tile', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_c', 
                         lvl2_id = 'nuts3', 
                         data = df_uk_socdist_scaled,
                         ctrls = 'prev')

extract_results(models_c_socdist)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

$random_slope_ctrl_prev

$interaction_ctrl_main_prev

$interaction_ctrl_int_prev

compare_models(models_c_socdist)

NA
NA

social distancing ~ extraversion

models_e_socdist <-run_models(y = 'socdist_single_tile', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_e', 
                         lvl2_id = 'nuts3', 
                         data = df_uk_socdist_scaled,
                         ctrls = 'prev')

extract_results(models_e_socdist)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

$random_slope_ctrl_prev

$interaction_ctrl_main_prev

$interaction_ctrl_int_prev

compare_models(models_e_socdist)

NA
NA

social distancing ~ agreeableness

models_a_socdist <-run_models(y = 'socdist_single_tile', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_a', 
                         lvl2_id = 'nuts3', 
                         data = df_uk_socdist_scaled,
                         ctrls = 'prev')

extract_results(models_a_socdist)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

$random_slope_ctrl_prev

$interaction_ctrl_main_prev

$interaction_ctrl_int_prev

compare_models(models_a_socdist)

NA
NA

social distancing ~ neuroticism

models_n_socdist <-run_models(y = 'socdist_single_tile', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_n', 
                         lvl2_id = 'nuts3', 
                         data = df_uk_socdist_scaled,
                         ctrls = 'prev')

extract_results(models_n_socdist)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

$random_slope_ctrl_prev

$interaction_ctrl_main_prev

$interaction_ctrl_int_prev

compare_models(models_n_socdist)

NA
NA

Explore quadratic trends

prevalence ~ openness

models_o_covid_exp <-run_models(y = 'rate_day',
                         lvl1_x = 'time',
                         lvl2_x = 'pers_o',
                         lvl2_id = 'ut_area',
                         data = df_uk_scaled,
                         ctrls = 'exp')

extract_results(models_o_covid_exp)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_exp

$interaction_exp

compare_models(models_o_covid_exp)

NA

prevalence ~ conscientiousness

models_c_covid_exp <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_c', 
                         lvl2_id = 'ut_area', 
                         data = df_uk_scaled,
                         ctrls = 'exp')

extract_results(models_c_covid_exp)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_exp

$interaction_exp

compare_models(models_c_covid_exp)

NA

prevalence ~ extraversion

models_e_covid_exp <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_e', 
                         lvl2_id = 'ut_area', 
                         data = df_uk_scaled,
                         ctrls = 'exp')

extract_results(models_e_covid_exp)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_exp

$interaction_exp

compare_models(models_e_covid_exp)

NA

prevalence ~ agreeableness

models_a_covid_exp <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_a', 
                         lvl2_id = 'ut_area', 
                         data = df_uk_scaled,
                         ctrls = 'exp')

extract_results(models_a_covid_exp)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_exp

$interaction_exp

compare_models(models_a_covid_exp)

NA

prevalence ~ neuroticism

models_n_covid_exp <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_n', 
                         lvl2_id = 'ut_area', 
                         data = df_uk_scaled,
                         ctrls = 'exp')

extract_results(models_n_covid_exp)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_exp

$interaction_exp

compare_models(models_n_covid_exp)

NA

Create overview table

Define function to create overview tables

summary_table <- function(models, dv_name, prev=F){

  temp_df_ctrl_main <- NULL
  temp_df_ctrl_int <- NULL
  temp_df_ctrl_int_prev <- NULL
  
  for (i in models){
    results <- i %>% extract_results()
    
    results_ctrl_main <- results$interaction_ctrl_main_dem['lvl1_x:lvl2_x',]
    temp_df_ctrl_main <- temp_df_ctrl_main %>% rbind(results_ctrl_main)
    
    results_ctrl_int <- results$interaction_ctrl_int_dem['lvl1_x:lvl2_x',]
    temp_df_ctrl_int <- temp_df_ctrl_int %>% rbind(results_ctrl_int)
    
    if(prev){
      results_ctrl_int_prev <- results$interaction_ctrl_int_prev['lvl1_x:lvl2_x',]
      temp_df_ctrl_int_prev <- temp_df_ctrl_int_prev %>% rbind(results_ctrl_int_prev)
    }
        
  }
  
  names_ctrl_main <- paste0(dv_name, '~', c('o', 'c', 'e', 'a', 'n'), '*time', '_crtl_popdens')
  rownames(temp_df_ctrl_main) <- names_ctrl_main

  names_ctrl_int <- paste0(dv_name, '~', c('o', 'c', 'e', 'a', 'n'), '*time', '_crtl_popdens*time')
  rownames(temp_df_ctrl_int) <- names_ctrl_int

  if(prev){
    names_ctrl_int_prev <- paste0(dv_name, '~', c('o', 'c', 'e', 'a', 'n'), '*time', '_crtl_popdens*time_prev')
    rownames(temp_df_ctrl_int_prev) <- names_ctrl_int_prev
    
    sum_tab <- rbind(temp_df_ctrl_main, temp_df_ctrl_int, temp_df_ctrl_int_prev) %>% round(4)
  }else{
    sum_tab <- rbind(temp_df_ctrl_main, temp_df_ctrl_int) %>% round(4)
  }


  
  return(sum_tab)

} 

Create overview tables

# prevalence
models_prev <- list(models_o_covid, 
                    models_c_covid, 
                    models_e_covid, 
                    models_a_covid, 
                    models_n_covid)

sum_tab_prev <- summary_table(models_prev, dv_name = 'prev')

write.table(sum_tab_prev, quote=F)

Value Std.Error DF t-value p-value
prev~o*time_crtl_popdens 0.0448 0.0025 3127 17.8777 0
prev~c*time_crtl_popdens -0.0386 0.0025 3127 -15.1522 0
prev~e*time_crtl_popdens 0.0372 0.0026 3127 14.5749 0
prev~a*time_crtl_popdens -0.0425 0.0025 3127 -16.8501 0
prev~n*time_crtl_popdens -0.0031 0.0026 3127 -1.1684 0.2428
prev~o*time_crtl_popdens*time 1e-04 0.0041 3126 0.0267 0.9787
prev~c*time_crtl_popdens*time 0.0039 0.0035 3126 1.1107 0.2668
prev~e*time_crtl_popdens*time 0.0047 0.0031 3126 1.5126 0.1305
prev~a*time_crtl_popdens*time -0.0025 0.0036 3126 -0.6789 0.4972
prev~n*time_crtl_popdens*time -0.0055 0.0024 3126 -2.2661 0.0235

# social distancing
models_socdist <- list(models_o_socdist, 
                       models_c_socdist, 
                       models_e_socdist, 
                       models_a_socdist, 
                       models_n_socdist)

sum_tab_socdist <- summary_table(models_socdist, dv_name = 'socdist', prev=T)

write.table(sum_tab_socdist, quote=F)

Value Std.Error DF t-value p-value
socdist~o*time_crtl_popdens 0.0068 0.0018 2791 3.7749 2e-04
socdist~c*time_crtl_popdens -0.006 0.0018 2791 -3.3046 0.001
socdist~e*time_crtl_popdens 0.0073 0.0018 2791 4.0335 1e-04
socdist~a*time_crtl_popdens -0.0071 0.0018 2791 -3.9228 1e-04
socdist~n*time_crtl_popdens -9e-04 0.0018 2791 -0.5053 0.6134
socdist~o*time_crtl_popdens*time -0.0012 0.003 2790 -0.4016 0.688
socdist~c*time_crtl_popdens*time 0.0015 0.0026 2790 0.548 0.5837
socdist~e*time_crtl_popdens*time 0.0026 0.0023 2790 1.1082 0.2679
socdist~a*time_crtl_popdens*time -8e-04 0.0027 2790 -0.2856 0.7752
socdist~n*time_crtl_popdens*time -0.0012 0.0018 2790 -0.6676 0.5044
socdist~o*time_crtl_popdens*time_prev -0.0016 0.0029 2789 -0.5482 0.5836
socdist~c*time_crtl_popdens*time_prev 6e-04 0.0025 2789 0.2296 0.8184
socdist~e*time_crtl_popdens*time_prev 0.0018 0.0022 2789 0.806 0.4203
socdist~a*time_crtl_popdens*time_prev -5e-04 0.0026 2789 -0.1883 0.8507
socdist~n*time_crtl_popdens*time_prev -6e-04 0.0017 2789 -0.3241 0.7459

Conditional random forest analysis

Extract slopes prevalence

# slope prevalence
df_uk_slope_prev <- df_uk %>% split(.$ut_area) %>% 
  map(~ lm(rate_day ~ time, data = .)) %>%
  map(coef) %>% 
  map_dbl('time') %>% 
  as.data.frame() %>% 
  rownames_to_column('ut_area') %>% 
  rename(slope_prev = '.')

# merge with control variables 
df_uk_slope_prev <- df_uk %>% select(-time, -rate_day) %>%
  distinct() %>% 
  inner_join(df_uk_slope_prev, by = 'ut_area') %>%
  drop_na()

head(df_uk_slope_prev)

NA

Extract slopes social distancing

# slope socdist
df_uk_slope_socdist <- df_uk_socdist %>% split(.$nuts3) %>%
  map(~ lm(socdist_single_tile ~ time, data = .)) %>%
  map(coef) %>%
  map_dbl('time') %>%
  as.data.frame() %>%
  rownames_to_column('nuts3') %>%
  rename(slope_socdist = '.')

# merge with control variables 
df_uk_slope_socdist <- df_uk_socdist %>% select(-time, -date, -socdist_tiles, -socdist_single_tile) %>%
  distinct() %>%
  inner_join(df_uk_slope_socdist, by = 'nuts3') %>%
  drop_na()

head(df_uk_slope_socdist)

NA

Explore distribution of slopes

df_uk_slope_prev %>% ggplot(aes(slope_prev)) + geom_histogram(bins = 100)

df_uk_slope_socdist %>% ggplot(aes(slope_socdist)) + geom_histogram(bins = 100)

df_uk_slope_prev

CRF prevalence ~ openness

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_o_fit_prev <- cforest(slope_prev ~ pers_o + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_prev[-1], 
                         controls = ctrls)

crf_o_varimp_prev <- varimp(crf_o_fit_prev, nperm = 1)
crf_o_varimp_cond_prev <- varimp(crf_o_fit_prev, conditional = T, nperm = 1)

crf_o_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_o_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

CRF prevalence ~ conscientiousness

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_c_fit_prev <- cforest(slope_prev ~ pers_c + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_prev[-1], 
                         controls = ctrls)

crf_c_varimp_prev <- varimp(crf_c_fit_prev, nperm = 1)
crf_c_varimp_cond_prev <- varimp(crf_c_fit_prev, conditional = T, nperm = 1)

crf_c_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_c_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF prevalence ~ extraversion

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_e_fit_prev <- cforest(slope_prev ~ pers_e + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_prev[-1], 
                         controls = ctrls)

crf_e_varimp_prev <- varimp(crf_e_fit_prev, nperm = 1)
crf_e_varimp_cond_prev <- varimp(crf_e_fit_prev, conditional = T, nperm = 1)

crf_e_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_e_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF prevalence ~ agreeableness

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_a_fit_prev <- cforest(slope_prev ~ pers_a + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_prev[-1], 
                         controls = ctrls)

crf_a_varimp_prev <- varimp(crf_a_fit_prev, nperm = 1)
crf_a_varimp_cond_prev <- varimp(crf_a_fit_prev, conditional = T, nperm = 1)

crf_a_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_a_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF prevalence ~ neuroticism

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_n_fit_prev <- cforest(slope_prev ~ pers_n + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_prev[-1], 
                         controls = ctrls)

crf_n_varimp_prev <- varimp(crf_n_fit_prev, nperm = 1)
crf_n_varimp_cond_prev <- varimp(crf_n_fit_prev, conditional = T, nperm = 1)

crf_n_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_n_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF social distancing ~ openness

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_o_fit_socdist <- cforest(slope_socdist ~ pers_o + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_socdist[-1], 
                         controls = ctrls)

crf_o_varimp_socdist <- varimp(crf_o_fit_socdist, nperm = 1)
crf_o_varimp_cond_socdist <- varimp(crf_o_fit_socdist, conditional = T, nperm = 1)

crf_o_varimp_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_o_varimp_cond_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF social distancing ~ conscientiousness

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_c_fit_socdist <- cforest(slope_socdist ~ pers_c + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_socdist[-1], 
                         controls = ctrls)

crf_c_varimp_socdist <- varimp(crf_c_fit_socdist, nperm = 1)
crf_c_varimp_cond_socdist <- varimp(crf_c_fit_socdist, conditional = T, nperm = 1)

crf_c_varimp_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_c_varimp_cond_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF social distancing ~ extraversion

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_e_fit_socdist <- cforest(slope_socdist ~ pers_e + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_socdist[-1], 
                         controls = ctrls)

crf_e_varimp_socdist <- varimp(crf_e_fit_socdist, nperm = 1)
crf_e_varimp_cond_socdist <- varimp(crf_e_fit_socdist, conditional = T, nperm = 1)

crf_e_varimp_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_e_varimp_cond_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF social distancing ~ agreeableness

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_a_fit_socdist <- cforest(slope_socdist ~ pers_a + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_socdist[-1], 
                         controls = ctrls)

crf_a_varimp_socdist <- varimp(crf_a_fit_socdist, nperm = 1)
crf_a_varimp_cond_socdist <- varimp(crf_a_fit_socdist, conditional = T, nperm = 1)

crf_a_varimp_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_a_varimp_cond_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF social distancing ~ neuroticism

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_n_fit_socdist <- cforest(slope_socdist ~ pers_n + airport_dist + males +
                          popdens + manufacturing + tourism +
                          health + academic + medinc + medage + conservative, 
                         df_uk_slope_socdist[-1], 
                         controls = ctrls)

crf_n_varimp_socdist <- varimp(crf_n_fit_socdist, nperm = 1)
crf_n_varimp_cond_socdist <- varimp(crf_n_fit_socdist, conditional = T, nperm = 1)

crf_n_varimp_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_n_varimp_cond_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))